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ABSTRACT 

This  report  addresses  the  development  of  a  model  for  blast  propagation  in  com¬ 
partmented  structures  with  progressively  failing  thin  structural  elements  like 
windows.  The  mathematical  model  is  based  on  the  assumption  that  the  frangi¬ 
ble  structural  features  (a)  are  thin  compared  to  the  characteristic  length  scale, 

(b)  do  not  deform  much  before  breakage,  and  (c)  fail  instantaneously  when  the 
corresponding  accumulated  damage  reaches  some  threshold.  The  presented 
numerical  model  employs  a  variation  of  the  Godunov  scheme  for  blast  propa¬ 
gation,  coupled  with  a  simple  governing  equation  for  the  accumulated  damage 
of  finite  elements  constituting  a  structural  feature  with  empirical  coefficients 
estimated  from  the  pressure-impulse  diagram.  When  the  accumulated  damage 
exceeds  some  threshold  value,  a  group  of  finite  elements  representing  a  desig¬ 
nated  entity  like  a  single  window,  is  forced  to  fail  thus  instantaneously  changing 
the  computational  domain  topology.  The  illustrative  numerical  simulations  in 
three  dimensions  demonstrate  that  the  model  behaves  reasonably  well  and  is 
capable  of  providing  rough  estimates  of  progressive  failure  and  blast  transfer 
within  and  between  compartmented  structures. 
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An  Accumulated  Damage  Model  for  Blast  Propagation  in 
Compartmented  Structures  with  Progressively  Failing  Thin 

Bulkheads 

Executive  Summary 

This  report  addresses  the  development  of  a  model  for  blast  propagation  in  compart¬ 
mented  structures  with  progressively  failing  thin  structural  elements  like  windows.  The 
mathematical  model  is  based  on  the  assumption  that  the  frangible  structural  features 
(a)  are  thin  compared  to  the  characteristic  length  scale,  (b)  do  not  deform  much  be¬ 
fore  breakage,  and  (c)  fail  instantaneously  when  the  corresponding  accumulated  damage 
reaches  some  threshold.  The  presented  numerical  model  employs  a  variation  of  the  Go¬ 
dunov  scheme  for  blast  propagation,  coupled  with  a  simple  governing  equation  for  the 
accumulated  damage  of  finite  elements  constituting  a  structural  feature  with  empirical 
coefficients  estimated  from  the  pressure-impulse  diagram.  When  the  accumulated  damage 
exceeds  some  threshold  value,  a  group  of  finite  elements  representing  a  designated  entity 
like  a  single  window,  is  forced  to  fail  thus  instantaneously  changing  the  computational 
domain  topology.  The  illustrative  numerical  simulations  in  three  dimensions  demonstrate 
that  the  model  behaves  reasonably  well  and  is  capable  of  providing  rough  estimates  of 
progressive  failure  and  blast  transfer  within  and  between  compartmented  structures  and 
hence  provides  useful  input  to  blast  effects  vulnerability  Sz  lethality  assessment  codes. 

The  concept  of  accumulated  damage  can  be  equally  applied  to  vulnerable  components 
placed  inside  a  structure.  The  component’s  accumulated  damage  normalised  by  its  critical 
damage  can  be  regarded  as  the  probability  of  failure  of  the  component.  This  approach 
replaces  the  notion  of  failure  probability,  which  is  very  hard  to  measure  and  thus  difficult 
to  validate  the  corresponding  model,  with  a  conceptually  simpler  notion  of  accumulated 
damage  governed  by  a  dynamic  equation. 

The  work  has  been  performed  under  the  National  Security  task  NS  07/002. 
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1  Introduction 

Interaction  of  blast  waves  with  structures  attracts  the  interest  of  many  scientists  and 
engineers  [Bangash  1993,  Jones  &  Brebbia  2006].  In  many  circumstances  the  most  impor¬ 
tant  question  concerns  the  estimation  of  blast  loads  on  the  structure  whose  structural  anal¬ 
ysis  is  done  after  solving  the  gas  dynamics  problem  in  the  known  flow  domain.  This  analy¬ 
sis  is  particularly  important  for  structural  engineers  dealing  with  protecting  critical  infras¬ 
tructure  from  the  explosive  events  [Rose  Sz  Smith  2002,  Smith  &  Rose  2002,  Zhou,  Hao  & 
Deeks  2005,  Remennikov  &  Rose  2005,  Smith  &  Rose  2006,  Remennikov  &  Rose  2007,  Clut- 
tera,  Mathis  &  Stahl  2007]. 

In  other  practically  important  problems  it  is  necessary  to  account  for  significant  defor¬ 
mation  of  structures  interacting  with  blast  waves.  In  this  case  the  governing  equations  of 
gas  dynamics  are  coupled  with  those  of  solid  mechanics  through  dynamically  changing  flow 
domain  that  becomes  part  of  solution  [Dowell  &  Hall  2001,  Gong  &  Andreopoulos  2009]. 
When  internal  structural  elements  are  fragile  and  have  relatively  little  strength  and  mass, 
such  as  conventional  glazing  elements,  the  effects  of  their  deformation  before  breakage 
and  of  the  debris  on  the  flow  field  can  be  neglected.  In  other  words,  when  pressure  load 
due  to  blast  damages  internal  bulkheads  and  windows,  the  solution  domain  changes  in¬ 
stantaneously  and  shock  waves  start  propagating  into  adjacent  compartments,  potentially 
causing  damage  to  other  structural  elements.  The  delay  of  breakage  of  internal  structural 
elements  results  in  a  more  complex  flow  pattern  due  to  reflected  shock  waves.  The  flow 
pattern  appears  to  be  very  sensitive  to  the  time  of  failure  which  is  quite  difficult  to  predict 
without  experiments. 

The  standard  approach  for  estimation  of  the  lifetime  of  structural  elements  is  based  on 
the  pressure- impulse  (P-I)  diagram  commonly  used  in  the  preliminary  design  of  protective 
structures  to  establish  safe  response  limits  for  given  blast-loading  scenarios  [Baker  1973, 
Baker  et  al.  1983,  Lloyd  1998,  Gelfand  &  Silnikov  2004,  Fallah  &  Louca  2007].  It  is 
assumed  that  the  time  history  of  the  loading  overpressure  p(t)  at  a  particular  element  can 
be  represented  by  a  typical  wave  profile  with  a  single  peak  overpressure  P  and  impulse  I 
defined  by 


P  =  max  pit) , 

(1) 

ta<t<  oc 

OO 

I  =  /  max  [pit),  0]  d t  , 

(2) 

ta 


where  t  is  time  and  ta  is  the  arrival  time  of  the  shock  wave.  Many  theoretical  studies 
employ  a  single-degree-of-freedom  model  with  the  pressure  time  history  approximated  by 
the  analytic  expression 

p(t)  =  P  exp  P  — —  ,  t  >  ta  .  (3) 

The  behaviour  of  a  structural  element  is  modelled  by  an  oscillating  mass  on  a  spring  driven 
by  the  time-dependent  force  (3),  and  the  maximum  deflection  of  the  mass  from  equilibrium 
can  be  explicitly  obtained  in  terms  of  P  and  /  [Baker  et  al.  1983].  Identifying  the  damage 
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level  with  the  maximum  deflection,  it  is  straightforward  to  plot  the  corresponding  iso¬ 
damage  curve  in  the  (P,  I)-plane,  which  is  called  the  P-I  diagram.  Typically,  the  P-I 
diagram  looks  like  a  shifted  hyperbola  approaching  the  pressure  asymptote  P  — >  Poo  as 
I  oo  (impulsive  loading  regime)  and  the  impulse  asymptote  I  — >  1^  as  P  — >  oo  (quasi¬ 
static  loading  regime). 

However,  the  P-I  diagram  concept  does  not  work  well  for  the  damage  prediction  of 
structural  elements  exposed  to  blast  waves  with  multiple  pulses,  since  the  pressure  time 
history  cannot  be  described  by  just  two  parameters,  P  and  I.  In  order  to  avoid  the 
complicated  problem  of  ambiguous  extracting  the  peak  overpressure  and  impulse  from  the 
pressure  time  history,  a  simple  model  based  on  the  concept  of  accumulated  damage  is 
proposed.  The  accumulated  damage  is  assumed  to  completely  define  the  internal  state 
of  the  structural  element,  and  therefore  its  rate  of  change  in  time  becomes  a  function 
of  the  damage  itself  and  of  the  immediate  environment  modelled  by  the  gas  states  at 
both  sides  of  the  thin  structural  element  [Antanovskii  2008  a,  Antanovskii  2009].  This 
function  has  to  be  phenomenologically  defined  or  experimentally  measured  after  providing 
a  rigorous  definition  of  the  accumulated  damage,  for  example,  by  identifying  it  with  the 
concentration  of  micro-fractures.  The  critical  damage  at  which  the  element  or  a  group 
of  elements  breaks  has  to  be  experimentally  measured.  A  simple  analytical  expression 
for  the  rate  of  accumulated  damage  that  involves  a  single  empirical  parameter  called  the 
cut-off  overpressure  is  suggested.  This  approach  allows  one  to  estimate  both  the  cut-off 
overpressure  and  the  critical  damage  from  the  pressure  and  impulse  asymptotes  of  the 
corresponding  pressure-impulse  diagram.  Note  that  the  accumulated  damage  is  implicitly 
defined  from  the  postulated  expression  for  its  rate. 

A  numerical  scheme  based  on  a  Godunov- type  solver  [Godunov  1959]  that  employs 
the  exact  Riemann  solver  for  flux  calculation  [Toro  1999],  coupled  with  forward-time  inte¬ 
gration  of  the  damage  accumulation  equation,  is  described.  A  pseudo-validation  test  and 
some  illustrative  numerical  simulations  in  three  dimensions  are  conducted,  which  demon¬ 
strate  that  the  model  is  capable  of  capturing  the  essential  physics  and  can  be  used  in 
vulnerability  studies,  particularly  when  an  accurate  result  is  hard  to  predict  due  to  a  high 
level  of  uncertainties. 


2  Mathematical  model 


The  mathematical  model  consists  of  gas  dynamics  equations  governing  blast  propa¬ 
gation  inside  and  around  a  structure  divided  into  compartments  by  internal  rigid  walls 
and  frangible  interfaces  representing  thin  structural  features  like  glass  windows  and  bulk¬ 
heads.  The  thin  structural  elements  respond  to  blast  loading  at  discrete  times  when  the 
accumulated  damage,  governed  by  a  certain  damage  accumulation  equation,  exceeds  some 
threshold  value.  In  this  case  the  structural  element  is  forced  to  fail  thus  dynamically 
changing  the  topology  of  the  computational  domain. 

For  the  sake  of  simplicity,  air  and  detonation  products  are  modelled  by  a  polytropic 
gas,  and  the  effect  of  gravity  is  ignored.  A  more  sophisticated  model  has  to  consider  a 
multi-phase  flow  of  a  variable-composition  mixture  of  gases.  Bearing  in  mind  the  numerical 
finite-volume  method  to  be  employed,  it  is  worth  starting  with  the  integral  form  of  gas 
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dynamics  equations. 

Let  t  be  time,  p  gas  density,  v  velocity,  p  pressure,  and  e  specific  internal  energy. 
The  classical  conservation  laws  of  mass,  momentum  and  energy  for  gas  flow,  written  for 
arbitrary  control  volume  u  with  piecewise  smooth  boundary  du j  oriented  by  inward  normal 
unit  vector  n,  are  as  follows 

^-JpdV  =  J  pvndA,  (4) 

dJ  did 


—  /  p  v  dV  =  /  (pvv-n  +  pn)  d  A  , 


did 


dt  r,e^w 


dK  = 


did 


p  \  e  +  \  lvl2 )  +p 


v  •  n 


dA. 


(5) 

(6) 


Here  dV  and  d A  denote  the  volume  and  area  elements,  respectively. 


The  integral  conservation  laws  imply  the  following  Euler  equations  of  gas  dynamics 
written  in  the  conservative  form 


!  +  V.(pv)=0, 


d{pw) 


dt 


+  V-(pv(g)v+pG)  =  0, 


d_ 

dt 


1 


P  1  e  +  -|v| 


+  V- 


P  1  e+  2|v| 


+  p 


v  >  =  0  . 


(7) 

(8) 

(9) 


where  V  denotes  the  gradient  operator,  <8)  the  tensor  product,  and  G  the  metric  tensor. 
The  balance  laws  are  completed  with  the  constitutive  equation  for  ideal  polytropic  gas 
[Courant  &  Friedrichs  1948] 


P  =  (  7-  1  )pe 


(10) 


where  7  is  the  ratio  of  the  specific  heat  at  constant  pressure,  cp.  to  the  specific  heat  at 
constant  volume,  cv. 

Equations  (7)— (10)  constitute  a  self-contained  model  for  gas  flow.  The  model  requires 
imposing  initial  and  boundary  conditions  for  gas  density,  internal  energy  and  velocity. 

The  initial  conditions  for  blast  propagation  are  defined  using  the  hot-gas  balloon  model 
for  the  products  of  detonation.  For  example,  knowing  the  total  mass  and  chemical  energy 
of  an  explosive  charge,  it  is  straightforward  to  calculate  the  initial  density  and  specific 
internal  energy  in  an  initial  volume,  usually  a  spherical  balloon,  occupied  by  the  products 
of  detonation.  In  most  cases  it  suffices  to  set  mass  and  energy  densities  uniform  in  the 
balloon  and  assume  zero  initial  velocity.  Finally,  the  specific  heats  ratio  of  the  products  of 
detonation  is  approximated  by  that  of  the  air  (7  =  1.4)  to  avoid  the  complex  problem  of 
multi-phase  flow  whose  numerical  solution  exhibits  spurious  density  profiles  [Kami  1994], 
It  is  conceivable  that  the  far-held  how  of  gas  is  mainly  inhuenced  by  the  total  mass  and 
energy  released  from  explosion. 
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Some  attention  should  be  paid  to  proper  modelling  of  boundary  conditions.  The  type 
of  boundary  conditions  depends  on  the  characteristics  of  the  hyperbolic  system  of  gas 
dynamics  equations.  In  our  situation  only  reflective  and  transmissive  boundary  conditions 
will  be  used,  which  are  numerically  modelled  by  introducing  a  layer  of  ghost  cells  just 
behind  the  boundary  of  the  computational  domain  with  an  adjusted  state  of  fictitious  gas 
[Oran  &  Boris  1987].  Clearly,  when  a  structural  element  fails,  the  boundary  conditions 
change  abruptly. 

The  interface  state  of  a  thin  structural  feature  is  completely  described  by  accumulated 
damage  5  considered  a  function  of  time  and  position  point  at  the  interface.  If  5  does 
not  exceed  some  threshold  value  5*,  the  interface  element  is  assumed  rigid  thus  reflecting 
shock  waves  from  both  sides  of  the  interface.  However,  when  5  >  6*,  the  interface  element 
is  forced  to  disappear  with  the  effect  of  changing  the  reflecting  boundary  condition  to 
natural  condition  between  the  adjacent  control  cells.  Since,  by  definition,  the  state  of  a 
system  determines  its  evolution,  the  rate  of  change  of  5  in  time  t  depends  on  5  itself  and 
the  immediate  state  of  environment.  This  results  in  the  following  ordinary  differential 
equation 

^  =  Q  (5, p+, p-,6+,6-)  (11) 

where  Q  is  some  non-negative  function,  and  p±  and  9±  denote  the  gas  pressure  and  absolute 
temperature  at  both  sides  of  the  structural  element,  respectively.  The  function  Q  can  be 
replaced  with  a  differential  operator  to  model  crack  propagation  along  the  structure. 

Assuming  the  pressure  difference  across  the  thin  structure  to  be  the  main  damage- 
causing  mechanism,  the  following  expression  for  damage  accumulation  is  employed 

Q  =  max(|p+  -p-\  -p*,  0)  (12) 

where  p*  is  some  cut-off  overpressure  depending  on  the  material  and  dimensions  of  the 
structure.  It  is  clear  that,  if  the  overpressure  \p+  —  p_  |  is  less  than  p*,  the  damage  rate  (12) 
vanishes  and  hence  no  increase  of  the  accumulated  damage  is  observed. 

To  the  first  approximation  the  parameters  p*  and  5*  can  be  estimated  from  the 
pressure- impulse  diagram  as  p*  =  P0 0  and  (5*  =  Ia 0.  Indeed,  6  and  5*  have  the  di¬ 
mension  of  impulse  defined  by  the  time  integral  of  pressure  according  to  Equation  (2).  It 
is  conceivable  that,  if  the  overpressure  |p+  —  p_|  does  not  exceed  the  value  of  the  pressure 
asymptote  P0 c,  there  is  no  structural  failure  irrespective  of  the  impulse.  Therefore,  the 
value  of  the  pressure  asymptote  has  the  meaning  of  the  cut-off  overpressure  p* .  Similarly, 
if  the  impulse  is  less  than  the  value  of  the  impulse  asymptote,  the  structural  element  does 
not  fail  irrespective  of  the  peak  overpressure.  For  high  pressure  compared  to  p*,  the  im¬ 
pulse  is  close  to  the  accumulated  damage  6  as  can  be  seen  from  Equations  (11)  and  (12), 
and  therefore  the  impulse  asymptote  can  be  identified  with  5*. 

It  is  worthwhile  noting  that  Equation  (11),  with  given  expression  for  the  right-hand 
side  (12),  implicitly  defines  the  accumulated  damage  5  by  assuming  that  initially  5  =  0. 
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3  Numerical  algorithm 


The  employed  numerical  model  involves  simultaneous  solution  of  the  Euler  equa¬ 
tions  (7)-(10)  and  damage  accumulation  equation  (11)  with  the  right-hand  side  (12). 
The  gas  dynamics  equations  are  solved  by  an  explicit  Godunov-type  conservative  scheme 
[Godunov  1959]  based  on  the  finite- volume  method  coupled  with  the  exact  Riemann  solver 
for  flux  evaluation  between  the  adjacent  control  cells  [Toro  1999].  As  a  by-product  of  the 
Riemann  solver,  the  gas  state  at  each  facet  between  two  control  cells  in  contact,  including 
the  boundary  and  ghost  cells,  is  calculated,  which  is  used  to  evaluate  the  damage  rate 
from  Equation  (12).  The  accumulated  damage  5  is  calculated  at  each  boundary  facet 
representing  a  structural  feature  using  the  explicit  forward-time  Euler  scheme.  When  the 
maximum  value  of  5  over  a  group  of  facets  representing  the  structural  feature  exceeds  <5*, 
all  the  facets  of  the  group  fail  thus  changing  the  internal  boundary  conditions. 

In  order  to  uniformly  handle  different  types  of  mesh,  either  structured  or  unstruc¬ 
tured,  the  control  cells  and  facets  separating  them  are  represented  by  a  directed  graph 
[Antanovskii  20086].  The  vertices  and  arcs  of  the  directed  graph  correspond  to  the  control 
cells  and  facets,  respectively.  The  graph  is  directed  because  each  facet  must  be  oriented. 
This  identification  is  possible  because  every  facet  has  exactly  two  adjacent  cells,  since 
every  cell  attached  to  the  boundary  has  a  corresponding  ghost  cell.  The  directed  graph 
changes  dynamically  whenever  any  facet  fails. 

Let  us  introduce  the  volumetric  density  of  momentum,  u  =  pv,  and  the  volumetric 
density  of  total  energy,  e  =  p  (e  +  \  |v|2j.  In  terms  of  these  conserved  variables  the 
conservation  laws  (4)-(6)  take  the  succinct  form 


l  J [p- u’ el  iV  =  / 

d)  did 


1  £+p 

u  n,  -  uu  n+pn,  - u  •  n 

P  P 


dA 


where 


(13) 


(14) 


according  to  the  equation  of  state  (10).  Let  {w*}  (i  =  1, . . . ,  Nc )  constitute  tessellation  of 
the  computational  domain  V  where  Nc  is  the  total  number  of  control  cells.  The  state  of 
the  gas  dynamics  flow  is  approximated  by  the  cell  average  values 

si  =  y  J  [p,  u,  e]  dE  (15) 

Wf 


and  the  fluxes  by  the  surface  integrals 


duii 


1  .  £+p 

u  n,  — uu  n+pn,  - u  •  n 

P  P 


dA 


where 


(16) 


Vi  = 


bJi 


(17) 
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is  the  cell  volume.  According  to  the  conservation  laws  (13)  and  the  equation  of  state  (14), 
the  following  system  of  ordinary  differential  equations  results 


d Si 
d  t 


=  Ri  ■ 


(18) 


In  order  to  complete  this  system  of  dynamic  equations  one  needs  to  express  the  rates  Rj 
in  terms  of  the  states  S ).  This  is  done  by  employing  the  Riemann  solver  which  allows 
one  to  explicitly  determine  the  fluxes  between  two  cells  in  contact,  provided  that  the  gas 
states  in  the  cells  are  approximated  by  constant  values  and  the  cells  are  separated  by  a 
planar  patch.  In  this  case  the  cell  average  values  St  are  approximated  by  the  first  order 
quadrature  formulae 


Ri  —  \Pii  Uj,  £j]  . 


(19) 


Let  us  assume  that  the  collection  {cj*}  (i  =  1, . . . ,  Nc)  is  composed  of  polyhedra,  poten¬ 
tially  of  different  type,  and  let  {pk]  (k  =  1, . . . ,  Nf)  be  the  collection  of  all  distinct  facets 
of  the  polyhedra.  The  ghost  cells  are  already  included  in  the  collection  of  polyhedra.  The 
state  of  the  ghost  cell  must  be  specified  to  model  the  appropriate  boundary  conditions 
[Oran  &  Boris  1987].  Denote  the  area  of  facet  ipk  by  Aj~,  select  a  unit  normal  vector  n^, 
and  introduce  a  IVy-by-2  facet-to-cell  connectivity  matrix  a  with  entries  defined  as  follows. 
If  facet  (fk  belongs  to  the  intersection  dcotl  n  du ;j2  and  points  from  c uq  to  then  set 
a(k,  1)  =  i\  and  a(k,  2)  =  %2 .  Mathematically,  the  cell  connectivity  matrix  a  defines  a 
directed  graph  (V,£)  where  vertices  V  are  cells  (including  the  ghost  cells)  and  edges  £  are 
facets. 

The  numerical  algorithm  is  straightforward.  First  update  the  states  of  the  ghost  cells 
according  to  the  imposed  boundary  conditions.  For  example,  for  an  ideally  reflective  wall 
or  symmetry  plane,  set  the  ghost  cell  state  equal  to  that  of  the  neighbouring  boundary 
cell  except  for  the  normal  component  of  momentum  which  must  have  the  opposite  sign. 
For  transmissive  boundary  condition,  set  the  ghost  cell  state  identical  to  the  boundary  cell 
state.  Then,  for  given  S)  at  some  instant  t,  calculate  the  cell  values  of  velocity  Vj  =  Uj/p* 
and  pressure  pi  from  Formula  (14).  Zero  out  the  array  of  rates  Ri.  For  each  facet  tpk  solve 
the  Riemann  problem  [Toro  1999]  using  the  given  gas  states  in  cells  un  and  Wj2  where 
i\  =  a(k,  1)  and  %2  =  a(k,  2).  Knowing  the  gas  state  at  the  facet  calculate  the  fluxes  of 
mass,  momentum  and  energy  along  the  unit  normal  vector  n^,  multiply  by  the  facet  area 
Ak,  and  add  to  the  rate  R*2  but  subtract  from  the  rate  R.lx  according  to  the  choice  of  the 
normal  n^.  Finally,  divide  the  calculated  rates  Ri  by  volumes  Vi.  This  procedure  defines 
the  cell  rates  {i?*}  as  a  function  of  the  cell  states  {S)}. 

The  accumulated  damage  is  approximated  by  a  constant  value  5k  defined  at  each  facet 
Pk  representing  the  thin  structural  feature.  Equations  (11)  and  (12)  naturally  transform 
to  the  system  of  equations 

^  =max(|p+-p-|  -p*,0)  (20) 

where  p ^  are  the  interface  pressures  at  both  sides  of  the  facet  obtained  from  the  solution 
of  the  Riemann  problems  as  a  by-product. 

As  a  result  of  this  procedure  one  ends  up  with  a  non-linear  system  of  ordinary  dif¬ 
ferential  equations  that  can  be  solved  by  the  explicit  Euler  scheme.  The  size  of  the  time 
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step  is  selected  according  to  the  Courant-Friedrichs-Lewy  (CFL)  condition  which  requires 
estimation  of  velocity  and  sound  speed  at  the  cell  interfaces.  Physically,  this  condition 
limits  the  time  step  to  the  value  small  enough  not  to  allow  Riemann  waves  to  propagate 
more  than  half  size  of  a  control  cell. 

Two  different  approaches  for  mesh  representation  are  considered.  In  the  first  approach 
the  mesh  is  represented  by  a  global  directed  graph  whose  vertices,  edges  and  the  connec¬ 
tivity  matrix  change  dynamically  whenever  any  facet  fails.  This  requires  to  implement  a 
quite  complicated  logic  as  the  control  cells  of  new  compartments  have  to  be  dynamically 
added  that  changes  the  numerical  model  size  and  may  cause  memory  problems.  In  the 
second  approach  the  mesh  is  represented  by  a  collection  of  directed  graphs,  one  directed 
graph  per  compartment.  The  allocated  memory  does  not  change  during  program  execu¬ 
tion  as  different  directed  graphs  exchange  information  through  their  own  ghost  cells.  The 
calculations  are  performed  only  for  those  compartments  which  are  affected  by  the  blast. 
Though  the  number  of  cells  is  slightly  increased,  the  second  approach  is  more  flexible  as 
the  coding  logic  is  considerably  simplified,  and  at  the  same  time  it  allows  one  to  imple¬ 
ment  a  parallel  multi-threaded  or  multi-CPU  code,  because  after  updating  the  states  of 
the  ghost  cells  of  the  directed  graphs  the  Riemann  solvers  can  be  executed  on  different 
processors  independently. 

The  described  numerical  algorithm  is  implemented  in  DBlast,  a  DSTO  software  for  sim¬ 
ulation  of  blast  loads  on  structures.  A  MATLAB®  graphical  interface  for  post-processing 
is  developed. 


4  Simulation  results 

All  the  simulations  are  conducted  for  air  (cv  =  716.46  J/kg/K,  cp  =  1003.51  J/kg/K) 
initially  kept  at  atmospheric  pressure  p  =  101.325  kPa  and  absolute  temperature  8  =  288  K 
(15°  C).  It  is  assumed  that  an  explosive  charge  of  trinitrotoluene  (TNT)  is  initiated  from 
its  centre,  and  all  the  chemical  reactions  are  completed  when  the  spherical  detonation  wave 
reaches  the  charge  boundary.  Therefore,  the  initial  density  of  the  products  of  detonation 
is  equal  to  the  density  of  the  condensed  TNT  explosive,  p  =  1.6  x  103kg/m3,  and  the 
initial  internal  energy  is  the  chemical  energy  released,  e  =  4.52  x  106  J/kg. 


4.1  Model  validation 

The  pure  gas  dynamics  solver  of  DBlast  has  been  successfully  validated  against  several 
benchmark  solutions,  such  as  Sod’s  shock  tube  problem  [Sod  1978],  along  the  lines  of  the 
gas  dynamics  solver  previously  developed  in  MATLAB®  [Antanovskii  20086].  In  addition 
to  this,  the  validation  of  DBlast  against  Air3D  developed  by  Rose  [2006]  and  ConWep 
(Conventional  Weapons  Effects  Calculator)  has  been  conducted.  A  spherically  symmetric 
blast  from  a  spherical  100  kg  charge  of  TNT  has  been  simulated  in  three-dimensions  by 
DBlast  and  Air3D  using  exactly  the  same  mesh  of  cubic  cells.  The  computational  domain  is 
a  cube  of  edge  length  5  m,  and  the  cell  size  is  25  mm.  The  charge  is  centred  at  one  corner  of 
the  cube,  the  origin  (0,  0,  0)  of  the  Cartesian  coordinate  system  (x,  y,  z)  aligned  with  three 
cube  edges,  and  the  three-dimensional  solution  has  three  planes  of  symmetry  x  =  0,  y  =  0, 
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and  z  =  0,  where  reflective  boundary  conditions  are  imposed.  Transmissive  boundary 
conditions  are  imposed  at  the  opposite  faces  of  the  cube.  The  pressure  gauge  points 
are  specified  at  the  space  diagonal  of  the  cube  at  4  m  and  5  m  stand-off  distances.  The 
pressure  time  history  shown  in  Figure  1  demonstrates  farely  good  agreement  with  Air3D 
and  ConWep.  It  is  worthwhile  noting  that  Air3D  uses  a  different  numerical  algorithm, 
whereas  ConWep  is  based  on  real  empirical  data. 

A  pseudo-validation  test  described  in  Example  4  (Testing  Glazing  in  a  Cubicle  Struc¬ 
ture)  of  Air  3D  users’  guide  [Rose  2006]  has  been  simulated  by  DBlast.  The  laminated  glass 
pane  of  dimensions  1.25  m  by  1.55  m  by  7.5  mm  has  been  selected  as  the  glazing  element. 
As  is  indicated  by  Johnson  [1999],  the  4.536  kg  charge  of  TNT  at  8  m  stand-off,  detonated 
1  nr  above  the  ground,  is  sufficient  to  cause  breakage  of  the  laminated  glass  window  in  a 
cubicle  structure  with  frontal  dimensions  2.15  nr  by  2.15  nr  and  length  3  nr,  whereas  3  kg,  at 
the  same  stand-off,  is  not.  The  simulation  results  obtained  with  the  use  of  Air3D  confirm 
this  assertion. 

Table  1:  P-I  data  for  a  laminated  glass  pane  of  dimensions  1.25  nr  by  1.55  nr  by  7.5nrnr. 


P  (kPa) 

2,000 

100 

70 

46 

30 

23 

20 

18 

17.5 

16.7 

16 

I  (kPa  X  ms) 

175 

180 

185 

200 

250 

300 

400 

500 

700 

1,000 

10,000 

The  corresponding  P-I  diagram  for  the  glazing  element  (taken  from  [Rose  2006])  is 
tabulated  in  Table  1  and  plotted  in  Figure  2.  It  is  clear  that  the  values  for  the  asymptotes 
of  the  P-I  diagram  are  as  follows 

p*  =  16kPa,  =  175  kPa  x  nrs  .  (21) 

The  geometry  of  the  computational  domain  is  shown  in  Figure  3,  where  the  solution  is 
assumed  to  have  two  planes  of  symmetry,  x  =  0  and  y  =  0.  The  obtained  results  simulated 
by  DBlast  have  predicted  the  same  effect  as  Air3D,  and  even  the  time  of  window  breakage 
is  virtually  identical.  The  pressure  time  history  at  the  centre  of  the  glass  window  and  the 
accumulated  damage  normalised  by  the  critical  damage  4*  are  depicted  in  Figure  4  for 
3  kg  charge  and  Figure  5  for  4.536  kg  charge. 


4.2  Blast  propagation  inside  a  building 

The  geometry  for  the  simulated  three-dimensional  structure  is  shown  in  Figure  6. 
The  two-room  building  contains  seven  windows,  one  of  which  is  internal,  of  the  same 
dimensions  as  indicated  in  Table  1,  so  the  cut-off  overpressure  and  critical  damage  are 
given  by  Formulae  (21).  The  walls  are  0.5  m  thick,  and  any  doors,  if  present,  are  shut.  A 
spherical  15  kg  charge  of  TNT  explosive  is  placed  at  position  x  =  4  m,  y  =  1  m,  z  =  0.5  m 
of  the  depicted  Cartesian  coordinate  system.  In  Figure  6  the  charge  can  be  seen  through 
two  windows  as  a  red  sphere.  The  computational  domain  of  this  example  consists  of 
7,  504,  530  cubical  cells  (excluding  ghost  cells)  and  22,  762,  238  square  faces.  It  suffices  to 
choose  the  CFL  number  equal  to  0.5  to  achieve  numerical  stability  that  resulted  in  3,  016 
cycles  to  complete  the  simulations  for  the  final  20  ms  physical  time. 

The  pressure  time  history  at  two  gauge  points  located  at  the  middle  of  each  room 
(floor  level)  is  shown  in  Figure  7.  Here  Room  1  refers  to  the  larger  room  where  the 
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charge  is  positioned,  whereas  Room  2  is  the  smaller  one.  The  pressure  contours  in  the 
vertical  (y  =  3  m)  and  horizontal  (z  =  1  nr)  cross-sections  are  displayed  in  Figures  8-15 
for  increasing  instances  of  time.  It  is  seen  that  the  windows  fail  in  turn,  thus  allowing  the 
blast  to  propagate  into  the  adjacent  room  (Room  2)  and  outside  of  the  building. 


5  Discussion 

The  described  model  is  significantly  simpler  as  compared  to  more  sophisticated  fluid- 
structure  interaction  models.  Nevertheless,  it  is  capable  of  capturing  the  effect  of  pro¬ 
gressive  failure  of  thin  structural  elements,  like  windows  and  light  bulkheads,  followed  by 
blast  transfer  into  adjacent  compartments.  In  many  circumstances  there  is  little  informa¬ 
tion  about  the  details  of  a  compartmented  structure  and  the  charge  characteristics,  and 
therefore  only  rough  estimates  for  the  level  of  damage  of  vulnerable  components  hidden 
inside  the  structure  is  required.  In  this  case  it  is  worthwhile  considering  as  simple  model 
as  possible  to  avoid  unnecessary  handling  of  uncertain  details. 

The  pressure-impulse  diagrams  (iso-damage  curves)  are  used  for  the  prediction  of  dam¬ 
age  of  other  materials  not  necessarily  glasses,  and  for  the  evaluation  of  the  response  of 
various  components  to  blast  loading.  However,  the  simplified  assumption  that  the  struc¬ 
tural  feature  completely  disappears  when  the  accumulated  damage  reaches  some  threshold 
is  applicable  to  sufficiently  thin  elements  composed  of  a  low-density  material. 

Note  that  the  proposed  model  does  not  take  into  account  the  dynamics  of  debris  and 
the  associated  exchange  of  momentum  and  energy  between  the  blast  and  fragments.  This 
shortcoming  can  be  overcome  by  considering  a  more  complex  model  for  the  structural 
element.  For  example,  it  can  be  modelled  by  a  layer  of  three-dimensional  rigid  elements, 
not  necessarily  of  zero  thickness,  until  it  breaks.  The  break  point  has  to  be  determined 
by  a  similar  criterion  for  accumulated  damage.  After  breaking,  the  material  of  the  rigid 
elements  is  replaced  with  a  fictitious  gas  of  high  density  and  of  different  equation  of  state 
(i.e.  different  specific  heats  ratio).  The  effect  of  debris  can  be  modelled  by  the  fictitious 
gas,  and  the  coupled  dynamics  of  blast  and  fragments  will  be  governed  by  a  multi-phase 
flow  model.  The  enhanced  model  does  not  track  each  individual  fragment  as  the  fictitious 
gas  represents  only  the  average  mass  of  fragments.  The  development  of  the  described 
multi-phase  model  for  blast  and  fragment  synergy  will  be  the  subject  of  a  separate  study. 

The  concept  of  accumulated  damage,  widely  used  in  fatigue  theory,  can  be  equally 
applied  to  vulnerable  components  placed  inside  a  structure.  The  component’s  accumulated 
damage  normalised  by  its  critical  damage  can  be  regarded  as  the  probability  of  failure  of  the 
component.  This  approach  replaces  the  notion  of  failure  probability,  which  is  very  hard 
to  measure  and  thus  difficult  to  validate  the  corresponding  model,  with  a  conceptually 
simpler  notion  of  accumulated  damage  governed  by  a  dynamic  equation. 
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Pressure  time  history  at  4  m  stand-off 


Pressure  time  history  at  5  m  stand-off 


Figure  1:  Validation  of  pure  gas  dynamics  solver  of  D  Blast  against  Air3D  and  ConWep. 
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Figure  2:  P-I  diagram  for  a  laminated  glass  pane  as  in  Table  1. 
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Glazing  element  in  cubicle  (3  kg  TNT  charge) 


x  (m)  y(m) 


Figure  3:  Computational  domain  for  testing  a  window  element  in  a  cubicle. 
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Gauge  point  at  the  middle  of  the  glazing  element 


Glazing  element 


Figure  f:  Pressure  time  history  at  the  centre  of  the  window  element  in  a  cubicle  and 
normalised  damage  for  3  kg  TNT  charge. 
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Gauge  point  at  the  middle  of  the  glazing  element 


Glazing  element 


Figure  5:  Pressure  time  history  at  the  centre  of  the  window  element  in  a  cubicle  and 
normalised  damage  for  4.536  kg  TNT  charge. 
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Blast  propagation  inside  a  building 


x  (m) 
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Figure  6:  Building  geometry. 
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Middle  of  Room  #1  (floor  level) 


Middle  of  Room  #2  (floor  level) 


Figure  1:  Pressure  time  history  at  the  middle  of  Room  1  and  Room  2  (floor  level). 
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Pressure  (kPa)  at  t  =  0.8  ms  (z  =  1  m) 
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Figure  8:  Pressure  contours  in  the  horizontal  and  vertical  cross-sections  at  t  =  0.8  ms. 
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Pressure  (kPa)  at  t  =  1 .4  ms  (z  =  1  m) 
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Figure  9:  Pressure  contours  in  the  horizontal  and  vertical  cross-sections  at  t  =  1.4  ms. 
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Figure  10:  Pressure  contours  in  the  horizontal  and  vertical  cross-sections  at  t  =  1.8  ms. 
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Figure  11:  Pressure  contours  in  the  horizontal  and  vertical  cross-sections  at  t  = 


2.1  ms. 
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Figure  12:  Pressure  contours  in  the  horizontal  and  vertical  cross-sections  at  t  =  4ms. 
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Figure  13:  Pressure  contours  in  the  horizontal  and  vertical  cross-sections  at  t  =  6  ms. 
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Figure  14 '■  Pressure  contours  in  the  horizontal  and  vertical  cross-sections  at  t  =  8  ms. 
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Figure  15:  Pressure  contours  in  the  horizontal  and  vertical  cross-sections  at  t  =  10  ms. 
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